clear;
doping=12.5;
t=0.6;
U=8;
rootPath = ['Hubbard_U8_doping',num2str(doping),'/'];
cptSpectrum = load([rootPath,'AkwCPTSpectra.txt']);
enrgList = load([rootPath,'enrgList.txt']); 
specInfo = load([rootPath,'specInfo.txt']); mu = specInfo(4);

momList = 1:size(cptSpectrum,1);
if size(cptSpectrum,1) == 501
    momList = [1:200 200+sqrt(2)*(1:100) (200+sqrt(2)*100)+(1:201)];
end

%Fermi Dirac function
kb=8.6e-5;
FD=1./(exp((enrgList-mu)*t/(kb*20))+1);
GFD=normpdf(enrgList*t,median(enrgList*t),0.1);
FDG=median(diff(enrgList*t))*conv(GFD,FD,'same');

EDC0=cptSpectrum(1,:);
cptSpectrumG=cptSpectrum;

G=normpdf(enrgList*t,median(enrgList*t),0.15);
for i=1:length(momList)
    cptSpectrumG(i,:)=median(diff(enrgList*t))*conv(G,cptSpectrumG(i,:),'same').*FDG;
end
SpectrumGSym=[fliplr(cptSpectrumG'),cptSpectrumG(2:end,:)'];
figure('position',[100,100,400,400]);hold on;box on; 

S=SpectrumGSym.^0.7;

pcolor(0:200, (enrgList-mu)*t, S); box on; shading interp;
xlabel('{\itq}/a'); ylabel('Binding Energy (eV)');

load('MyColormaps_Terrain','mycmap');
terrain1=mycmap(1:350,:);
colormap(terrain1);

caxis([0.03 max(max(S))]);
set(gca,'tickdir','out'); 
ylim([-1.5 0.2]);
xlim([0 200]);
set(gca,'XTick',[0 100 200]);
set(gca,'XTickLabel',{'-\pi','0','\pi'});
set(gca,'PlotBoxAspectRatio',[1 1 2]); 
set(gca,'LineWidth',1);
set(gca,'FontSize',20);
% saveas(gcf,['CPT',num2str(round(-doping*2/64*100)),'.png']);
plot(0:200,0*(0:200),'--w');
% colorbar

